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-TITLE MTHSATANH 3; Single Precision Hyperbolic Arctangent 
"IDENT /2-003/ t Pile: MTHATANH MAR’ Edit: PDG2O0S © 


COPYRIGHT Wfe 1978, 1980, 1982, 1984 BY 
DIGITAL EQUIPMENT CORPORATION, MAYNARD, MASSACHUSETTS. 
ALL RIGHTS. RESERVED. 
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* THIS SOFTWARE 1S FURNISHED UNDER A LICENSE AND MAY BE USED AND COPLED 
* ONLY IN ACCORDANCE WITH THE TERMS OF SUCH LICENSE AND WITH THE 
* INCLUSION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER 
* COPIES THEREOF MAY NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY 
: OTHER PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE 1S” HEREBY 
* TRANSFERRED. 

* 

* THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE 
* AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT 
* CORPORATION. 
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DIGITAL ASSUMES NO RESPONSIBILITY FOR THE USE OR RELIABILITY OF ITS 
SOFTWARE ON EQUIPMENT WHICH IS NOT SUPPLIED BY DIGITAL. 
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FACILITY: MATH LIBRARY 
“ABSTRACT: 


MTHSATANH returns the single precte ioe hyperbolic arctangent of the 
single precision argument. The call is standard call-by-reference. 


VERSION: 2 
HISTORY: 
AUTHOR: 
Peter D Gilbert, 23-Jul-81: Version 2 


MODIFIED BY: 
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ORY ; 
8 47 -SBTTL HISTORY ; Detailed Current Edit History 
; VERSION 1 

000 


1-001 = Algorithm from PL/I math Library. 
Edit History for Version 02 of MTHSATANH 


QO Rewrite of PL/I version. July 1981 

1 = Change MOVZBL te CVTBL when accessing ty mer BAL OG. V. PDG 2-Dec-1981 
¢ = Change RSB to RET after error exit. PDG 

-005 - Repair problem with POLY instruction. PDG 1Gemar=1982 
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MTHSATANH 
2-003 


Se ate 


09 69 -SBTTL DECLARATIONS ; Declarative Part of Module 
62 ; 
és ; INCLUDE FILES: MTHJACKET.MAR 
9? 
308 ef ; EXTERNAL SYMBOLS: 
0 0 63 .DSABL GLOBAL 
00 0 «SHOW BINARY,CALLS,CONDITIONALS,DEF INITIONS, EXPANSIONS 
00 71 “EXTRN MTHS$K_{NVARGMAT 
ie ~EXTRN MTHSSSIGNAL 
09 7 -EXTRN MTHS$$AB_ALOG_V 
-EXTRN MTHSSAB_F_FH 
00 74 EXTRN MTHSS$AB-F FHT 
9000 fe 
0000 5 ; EQUATED SYMBOLS: 
000 78 ; 
0000 79 
0000 80 ; 
0000 81 ; MACROS: 
$83 
0000 Hn MACRO OPDEF xX, OP, SH 
0000 3s 85 OPDEF ADDX  *kOO@SH+OP,R'X,M'X 
0000 43-8 .OPDEF ADDX3 “*X01@SH+OP,R'X.R'X,W'X 
0000 87 OPDEF SUBX *x02@SH+OP,R'X.M'X 
0000 = s«&B OPDEF SUBX3 “XO3@SH+OP,R'X,R'X,W'X 
0000 836 89 -OPDEF MULX  “XO04@SH+0P,R'X,M'X 
0000 90 F MULX3 “XO5@SH+OP,R'X,R'X,W'X 
0000~=Ss(«91 OPDEF DIVX “X06aSH+OP,R'X,M'X 
0000 3¢ OPDEF DIVX3 “XO7@SH+OP,R'X.R'X,W'X 
0000 («9 OPDEF CVTWX  “XOD@SH+OP,RW,W'X 
0000 49-49 OPDEF POLYX  “X15@SH+OP.R'X,RW,AB 
0000 38 95 OPDEF Vx “x00D0,RL,WL ; MOV 
0000 «4% OPDEF MOVAX “XOODE.AL.WL > MOVAL 
$000 98 
000 48699 OPDEF FF, <*x0040>, 0 
000 .OPDEF ADDX  “*x06a0+*X0040,RF MF 
000 -OPDEF ADDX3 2X9180420060 ,RF RE WF 
000 OPDEF SUBX *x0200+*X0040,RF ,MF 
99 OPDEF SUBX3 “x03a0+*x 04 RE RE WE 
OPDEF MULX *x04a0+°X0040.RF .MF 
900 OPDEF MULX3 °X05a0+*X0040.RF.RF.WF 
000 OPDEF DIVX *x06a0+*X 040.RE ME 
909 .OPDEF DIVX3 *X07a0+*X0040_RF RF WE 
09 -OPDEF CVTWX  “xODa0+*x 040.RW. WF 
00 OPDEF POLYX  *X15a0+*xX0040.RF.RW,AB 
000 OPDEF MOVX “xOODO.RL,WL : MOVL 
000 OPDEF MOVAX  “XOODE,AL,WL : MOVAL 
200 100 
00000007 0 101 F_EXP = 7 ; Bit offset to exponent 
00000007 44 ! § X_EXP = 7 ; Bit offset to exponent 
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: PSECT DECLARATIONS: 


-PSECT _MTHSCODE PIC,SHR,LONG, EXE ,NOWRT 
; program section for math routines 


; OWN STORAGE: none 
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recision Hyperbolic Arc A ant 16- 
Ss; Declarative Part cf 


Modul 6=-S 
; CONSTANTS: 


ACMASK = “M<R2,R3,R4,R5> 


LN2_HI: .LONG 472003681 
LN2-LO: .LONG  “XBEBE333F 
LOGTAB1: 
.LONG “*XEABDBF2A 
;LONG  *XOCDD3F4D 
{LONG “XFFFOBF7F 
“LONG “XAAA73FAA 
.LONG *x0000C000 
: Remove this constant, 
rigsgsz. /LONG *x00000000 
LOGLEN = .-LOGTABI/4 = 
LOGTAB2: 
.LONG  *X6D943FCD 
“LONG  “*XAA91402A 
:LO “00004100 
LOGLEN2 = .-LOGTAB2/4 = 1 


AX/VMS 
MTHRTL. 


p=19B6 11:20:39 


v04-00 
H 


acro 
RCIJMTHATANH.MAR; 1 


; register cet mask and integer 
; overflow enable 


: (Hi 16 bits of ln2)*2**-7 
3; (Low bits of Ln2)*2*2-7 


(2). Generated using $9 eq. 
et. al. (sin(2a) = 1/732) 


AOOOOAN 
—MWSLvur O 


and do another multiply in-Line. 


; CO = .00000000 
; no. of floating point entries 


Constants for p(z*z). ag using 
eq. 6.35.11 of Hart et. al. (sin(2a) = 
(b - 1)/(b + 1) where b = 2**(1/7)) 

C2 = 0.40122664 

C1 = 0.66666514 

CO = 2.00000000 
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-SBTTL MTHSATANH - Single Precision Hyperbolic Arctangent 


FUNCTIONAL DESCRIPTION: 

ATANH = Single precision floating point function 

ATANH(X) is computed using the following approximation technique: 
If iXi >= 1.0, error. Otherwise 
Let (1#xX)/(1=X) = f ® (2e#n), where 1/2 <= f < 1 


If °. bs Ff eyed — or equal to 1 then 
stee =n- d Fl = 2ef, 


"eet N=n and F = f. 


If if = 1! < 2e*-5 then 
2tatanh(Xx) = NeLn(2) + W + WeP(W), 
mare Wo = (C14F)/C1<F)) e2een = 
and P is a polynomial of rine Feel, f= 5,D=9. 


* obtain FHI (roughly equel to F) from table look 

2tatanh(xX) = Ln€(14#x)/(1=x)) = Neln(2) + in(FHI) + Z*Q(2*Z), 
where Q is a pop yogaial of of afreree F=2,D=5, 
where Z = (F = HI) 
where F = (2e8-N) #199070) 


cy 


Z is senputer Oy: 

Z = D)/(C1=xX*D) 
AR, Y = FHI *2**N 
where D = (¥=-1)/(Y¥4#1) 

Note that Z may be computed in a variety of ways: 
Z = COV#K) = Y¥e1-X)I/0014xX) + ¥e(1=x)] 
Z=(1 +X -Y + XevyJ/C1 +X + Y¥ = X*V¥J 
Z={(1-Y¥ +X + Xey)/C1 + ¥ + X = X*Y¥J 
Z = CO1-¥) + MeCT4evI/OC1+¥) + X#(1-¥)) 


NOTE: The ayers ties |n(A=FHI) and Ln2 are used in the above 
oqust’ ons in two parts - a high part (containing the 
—— order bits) and a low part (containing the low 
er bits. In the code the h gh and _ pects of the 
constants are indicated by a _H suffix eqegees 
tivel The values were chosen a= thet N*LN2_HI + 
tn. FH “HI is exactly representable. 
CALLING SEQUENCE: 
atanh.wf.v = MTHSATANH(x.rf.r) 
INPUT PARAMETERS: 
X=4 ; Contents of x is the argument 
3 IMPLICIT INPUTS: none 


: OUTPUT PARAMETERS: 
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0 3 
: VALUE: Single precision hyperbolic arctangent of the argument 
5 : IMPLICIT OUTPUTS: none | 
0 0? ; COMPLETION CODES: none | 
0 03 : SIDE EFFECTS: | 
Be 1] ; Signals: MTHSK_INVARGMAT if !X! >= 1.0 with reserved operand in RO (copied to | 
0 \¢ 3 the signal mechanism vector CHF$L_MCH_RO/R1 by LIBS$SIGNAL). | 
88 13 ; Associated message is: ‘Floating Overfiow in math Library’’. Result is 
14; reserved operand -0.0 unless a user supplied (or any) error handler changes 
99 13 ; CHFSL_MCH_RO/R1. | 
00 13 : NOTE: This procedure diseoles floating point underflow and integer 
00 18 ; overflow, causes no floating overflow or other arithmetic traps, and 
B05 4 ; preserves enables across the call. 
0028 1 3 Note: This routine is written to avoid causing ony integer overflows, 
0028 2 ¢ : floating overflows, or floating underflows or divide by 0 conditions, 
0028 ; whether enabled or not. 
0028 4; 
0028 25 ;--- 
0028 6 
0134— 31 Ose f ERR: BRw ERROR 
003¢ boSe 9 -ENTRY MTHSATANH, ACMASK ; standard call-by-reference entry 
002D 230 ; disable DV (and FU), enable IV 
50 04 BC 00 one =) MOVX @X(AP), RO ; RO = arg 
52 08 50 43 0031 = SUBF3 RO, S*#1.0, R2 ; R2 = 1-x 
ee) aS ite 3 BLEQ 3; ATANH(X) is not defined for X>=1 
54 08 50 41 0037 5 ADDF3 RO, S*#1.0, R4 3: R4 = 14x | 
FB 13 0038 336 BLEQ. ERR ; ATANH(X) is not defined for X<=-1 
54 52 46 003D 37 DIVF2 R2, R4 ; R4 = approximation to (1+X)/(1-x) 
55 00000000'GF 9E 0040 38 MOVA G*MTH$$AB_ALOG_V, R5 
5 65 CO 0047 39 ADDL (R5), R5 ; RS = address of ALOG table 
53 54 OO7F 8F AB O04A 40 BICW #iar BxP R4, R3 ; R3 = Biased exponent 
53 4000 8F A p20 41 SUBW #°x4000, RS ; R3 = Unbiased exponent 
63 1 bee? +g BLEQ NEG_EXP ; Branch to processing for n=<0 
53 0080 8F Ag 0057 44 SUBW #iaF_EXP, R3 ; Exponent is positive, R3 =N=n- 1 
54 ¢? A Opec 45 SUBW 3, R4 7 RG =F = 
54 4 9A OO05F 46 MOVZBL R4, R4 ; R4 = index into ALOG table 
00000000 ane 47 IF NE, F_EXP=X_EXP 
06 48 Divwe #1a<F-EXP-X"EXP>, R3 ss: Shift R3 to scale X-floating 
7E 53 4D 006 50 CVTWX R3, -(SP) ; Push N onto the stack | 
55 6544 98 006 51 CVTBL © (R5)CR4), R5 ; RS = offset into FHI tables | 
3¢ 19 0069 3 SS. LN_1_PLUS_ : Branch to handle F close to 1 | 
55  00000000'GF4 DE oo8 5 MOVAX G*ATASSAB_F_FHICRSJ, R5 ; RS = Address of FHI 
54 85 00 Bf? 2 MOVX (R5)+, R4 ; R4 = FHI 
a076 2$ > Compute Z = (F = FHI)/(F + FHI) | 
th af 3 Z = tay - a ge (14#X) + Y¥e(1-XxX)] 
0076 58 ; Z=(1+X-=-¥ + Xev¥)/(1 + X + ¥ = Xe¥) 
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50 


55 00000000 ' GF 
Py) 


; Single Pre 
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54 53 9 
98 54 4 
4 8 40 
P) 4 46 
> 0 41 
P) Q 44 
5 08 40 
54 5246 
54 54 «45 
02 50 55 
50 54 4G 
CF 6€ 45 
52. 85 40 
50 652 0 =640 
45 
6 40 
50. 5200S 40 
0080 8F A2 
04 
54 53 A2 
54. 9A 
00000000 
7E 36530— 4D 
6544 98 
50.19 
45 DE 
65 00 
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Compute Z**2, 
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R2 


P(Z**#2) and Z2*P(7**2) 


(SP), LN2_LO, R2 
(R5)+, R2 
R2, RO 


(SP)+, LN2_HI, R2 
(RS), Re 


R2, RO 
#1ax_EXP, RO 


R3, R4 


F_EXP=X_EXP 
#1<FEXP=X"EXP>, R3 


R3, -(SP) 
(R5)CR4], RS 


(R5) 


(fF = yee + FHI) 


R4, RG, RO 
RO, #LOGLEN2, LOGTAB2 
R4, 0 


91:98:36 AX/VMS Mac 
) 


where Y = FHI*2**N, roughly equal 


ro V04-00 
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to (14#X)/(1=x) 


R4 = FH] * 2e*N = SFHI 
R2 = 1 = SFHI 

R4 = 1 + SFHI 

R2 = (1-SFHI)/(14+SFHI) = D 
R4 =D + X 

R2 =D * X 

R2 = 1 + Dex 

R4 = (D#X)/(14D*X) = 7 
RO = Z**2 

RO = P(Z*#2) 

RO = Z*P(7*#2) 


; Compute B = N*LN2_LO + LN_FHI_LO + 2#P(Z#Z) 


N*LN2_LO 
N*LN2_LO + LN_FHI_LO 


; Compute A = N*LN2_HI + LN_FHI_HI and ALOG(X) 


R2 = N*LN2_HI 

R2 = A = N*LN2_HI + LN_FHI_HI 
RO = A + B = ALOG(X) 

Divide by 2 

R4 2 Ff = 2f 

R4 = index into ALOG table 


Shift R3 to scale X-floating 


Push WN onto the stack 
RS = offset into FHI tables 


PLUS ; Branch to handle F close to 1 
AS$AB_F_FHICRSJ, RS ; RS = FHI 
» RG ; RG = FHI 


Address of 


z= #X) = YaCT=X)I/CC14+X) + ¥e(1-x)) 
Z = (C1 +X -Y¥ + xXevyjJ/(1 + X + ¥ = Xe¥) 
= FHI*2**N, roughly equal 


XY 
to (14#X)/(1-x) 


R4 = FHI * 2**N = SFHI 

R2 = 1 = SFHI 

R4 = 1 + SFHI 

R2 = (1-SFHI)/(1+SFHI) = D 
RG =D +X 

R2 = D* X 

R2 = 1 + Dex 

R4 = (D#X)/(14D*x) = 7 
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EE 1 | 
OFF 1? § Compute Z**2, P(Z**2) and 2*P(7**2) 
54 54 «645 E 18? MULX3 =R4, 3 RO = Z**#2 
FF24 CF 08 20 55 OOF 0 POLYX RO, ; HoctEN2, LOGTAB2 ; RO = P(Z*#2) | 
: 4 44 0088 1 MULX ; RO = Z2*P(Z*#2) | 
ore § | ; Compute B = N*LN2_LO + LN_FHI_LO + Z*P(2*Z) | 
52 FFOG CE 6 5 O0FB 5° MULX3 (SP), LN2_LO, R2 ; R2 = N*LN2_LO 
52 of 40 101 § ADDX = (R55 NS 3 RS = N*LN “LO + LN_FHI_LO 
50 «52)—40 194 ADDX  R2, RO : RO=8 
o107 ; ; Compute A = N*LN2_HI + LN_FHI_HI and ALOG(X) 
52 FEF4 CF oe 45 107 1 MULX3 (SP)+, LN2_HI, R2 3; R2 = N*LN2_HI 
$2 7? $6 010D 3 SUBX -(R5), R2 3; R2 = A = N®LN2_HI + LN_FHI_HI 
- . Soe g119 3 ADDX R2, R ; RO = A + B = ACOG(X) 
50 0080 8F A2 011 34 SUBW2 #1aXx_EXP, RO 3 Divide by 2 
04 0118 335 RET 
0119 336 
0119 tf ; 
0119 8 ; Special logic for F close to 1 
0119 39 ; 
0119 340 | 
0119 341 LN_1_PLUS: 
54 08 50 43 0119 325 SUBX3 RO, $*#1.0, R4 : = 1-x 
53 BS 011D 34 TSTw R3 : pet termine which way to calculate W | 
OF 13 O11F 344 Beat. 108 
54 10 54 47 0121 45 DIVX3 af See. 0. RG 3 R4 = 2/(1-x) 
54 98 42 b1S8 46 SUBX “#1, 0, 3 R4 = (14K) /01-x) 
6 § as 0128 47 SUBW 3: Scale R4 
54 08 #442 0128 348 SUBX m3, Ra R4 : RG = W 
09 «11 orgs 49 BRB 
54 50 54 47 0130 350 10$ DIVX3 R4, RO, RS 7 R4 = xf (1-x) 
54 0080 8F AO 0134 $2) ADDW #iax EXP, R4 ; R4 = = sex) = (14X)/(1-x) = 1 
FEC9 CF 04 54 55 0139 2¢ 20$: POLYX R4, FLOGLEN1, LOGTAB1 3; RO = aw 
50 54 44 O13F 5 MUL X R4, RO : Finish p congut tng Q(wW) 
50 54 44 10g 54 MULX R4, RO 3; RO = 
52 FEBA CFE 6€ 45 014 55 MULX3 (SP), LN2_LO, R2 : R2 = NOUN LO 
50 52 40 0148 356 ADD R2, RO : RO = N®*LN “0 + Wea(W) 
50 54 40 O14E 57 ADDX R4, RO ; RO = N®LN O + WeQ(W) + W 
6€ FEAB CF 44 0151 358 MULX LNO_HI, (SP) : (SP) = NeLNo_ 
50 ge 40 b128 59 ADDX (SP)+, R 3 RO = ALOG(X) 
50 0080 BF A2 015 60 SUBW2  #lax_EXP, RO : Divide by 2 | 
04 Bide o RET 
O15 68 | 3; X <= 0.0, signal error | 
00'S8F 9A O15F 65 ERROR: a #MTHSK_INVARGMAT, -(SP) ; condition value 
50 01 OF 78 0163 6 #15, #T, RO 3 RO = result = reserved operand -0.0 
0167 6 : goes © signal mechanism vector 
167 o8 : CHSt OncH -RO/R1) so error handler 
167 § 3 can modify the result. 
00000000'GF 01 FB 0167 0 CALLS #1, G*MTHSSSIGNAL 3; signal error and use real user's PC 
16E 71 ; independent of CALL vs JSB 
04 O16E 7Z RET 3; return - RO restored eo 
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Big | 32 END ; CHFSL_MCH_RO/R1 
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Symbol table 66-SEP-1984 11: 
ACMASK = £00008 f 
ERR 8 Rp 4 
ERROR 1SF R 1 
F_EXP = 00000007 
LN2_HI 000 8 fe 01 
LN2_LO siti 4 R 01 
LN_T_PLUS 0 Mtg ie 01 
LN_1 PLUS_W 44 8 01 
LOGLEN1 = 49444 04 
LOGLEN2 = 0000 it 
LOGTAB1 00000008 R 01 
LOGTAB2 0000001C R 01 
MTHSSAB_ALOG V  teeeeere =X = 00 
MTHSSAB_F FHT eeeeeeee =X 00 
MTHS$SIGNAL wererere = X 00 
MTHSATANH 00000028 RG 01 
MTHSK_INVARGMAT ‘*teeeeee x 00 
NEG_EXP QOOOOOBA R 01 
= 00000004 

X_EXP = 00000007 

tro w ee eee eo a wmmance a 

! Psect synopsis ! 

tomar wmowrocnr soe ncaa + 
PSECT name Allocation PSECT No. Attributes 
- «ABS, 00000000 ( 0.) 00 ¢ OQ.) NOPIC USR CON ABS LCL 
_MTHSCODE OOOOO16F (¢ 367.) O1¢ 1.) Pic USR CON REL LCL 

Denaniecinensenampescntbiintainabnineatettaeniaiits 
' Performance indicators ! 

Phase Page faults CPU Time Elapsed Time 
Initialization 33 00:00:00.08 00:00:00.72 
Command processing 126 00:00:00.67 00:00:04.14 
Pass 94 00:00:01.28 00:00:04.37 
Symbol table sort 0 00:00:00.01 00:00:00.01 
Pass 2 85 bs bs OF 00:00:03.66 
Symbol table output 3 BF 8 088 00:00:00.03 
Psect synopsis output 00:00:00.02 00:00:00.02 
Cross-reference output 4 $4 4 BFR 8-8 
Assembler run totals 34 0:00:03. 00:00:12.98 


The working set Limit was 900 pages. 

6708 bytes (14 pages) of virtual memory were used to buffer the intermediate code. 

There were 10 pages of symbol table space allocated to hold 20 non-local and 2 local symbols. 
434 source Lines were read in Pass 1, producing 11 object records in Pass 2. 

3 pages of virtual memory were used to define 2 macros. 


MTHRTL. 


NOSHR NOEXE NORD NOWRT NOVEC BYTE 
SHR EXE RD NOWRT NOVEC LONG 


VAX-11 Macro Run Statistics 


tewese mero mce roe reese ecmeosae} 


! Macro Library statistics ! 


Quncnes wernt enoeeeeea anna 


Macro Library name 
-$255$DUA28: CSYSLIBJSTARLET.MLB; 2 0 
0 GETS were required to define 0 macros. 


There were no errors, warnings or information messages. 
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cro V04-00 Page 18 
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MACRO/ENABLE=SUPPRESSION/DISABLE=(GLOBAL , TRACEBACK) /LIS=LIS$:MTHATANH/OBJ=O0BJ$:MTHATANH MSRC$:MTHJACKET/UPDATE=(ENHS:MTHJACKET) #MSRC 
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